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Abstract 

The  level  set  method  [31]  is  a  popular  technique  for  tracking  moving  interfaces 
in  several  disciplines  including  computer  vision  and  fluid  dynamics.  However,  de¬ 
spite  its  high  flexibility,  the  original  level  set  method  is  limited  by  two  important 
numerical  issues.  Firstly,  the  level  set  method  does  not  implicitly  preserve  the 
level  set  function  as  a  distance  function,  which  is  necessary  to  estimate  accurately 
geometric  features  s.a.  the  curvature  or  the  contour  normal.  Secondly,  the  level 
set  algorithm  is  slow  because  the  time  step  is  limited  by  the  standard  CFL  con¬ 
dition,  which  is  also  essential  to  the  numerical  stability  of  the  iterative  scheme. 
Recent  advances  with  graph  cut  methods  [4,  3]  and  continuous  convex  relaxation 
methods  [7,  5,  16]  provide  powerful  alternatives  to  the  level  set  method  for  image 
processing  problems  because  they  are  fast,  accurate  and  guaranteed  to  find  the 
global  minimizer  independently  to  the  initialization.  These  recent  techniques  use 
binary  functions  to  represent  the  contour  rather  than  distance  functions,  which  are 
usually  considered  for  the  level  set  method.  However,  the  binary  function  cannot 
provide  the  distance  information,  which  can  be  essential  for  some  applications  s.a. 
the  surface  reconstruction  problem  from  scattered  points  and  the  cortex  segmen¬ 
tation  problem  in  medical  imaging.  In  this  paper,  we  propose  a  fast  algorithm 
to  preserve  distance  functions  in  level  set  methods.  Our  algorithm  is  inspired  by 
recent  efficient  i1  optimization  techniques,  which  will  provide  an  efficient  and  easy 
to  implement  algorithm.  It  is  interesting  to  note  that  our  algorithm  is  not  limited 
by  the  CFL  condition  and  it  naturally  preserves  the  level  set  function  as  a  distance 
function  during  the  evolution,  which  avoids  the  classical  re-distancing  problem  in 
level  set  methods.  We  apply  the  proposed  algorithm  to  carry  out  image  segmen¬ 
tation,  where  our  methods  proves  to  be  5  to  6  times  faster  than  standard  distance 
preserving  level  set  techniques.  We  also  present  two  applications  where  preserving 
a  distance  function  is  essential.  Nonetheless,  our  method  stays  generic  and  can  be 
applied  to  any  level  set  methods  that  require  the  distance  information. 
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1  Introduction 


In  the  last  twenty  years,  the  level  set  method  (LSM)  of  Osher  and  Sethian  [31]  has 
become  a  popular  numerical  technique  for  tracking  moving  interfaces  in  computational 
geometry,  fluid  mechanics,  computer  graphics,  computer  vision  and  material  sciences. 
The  main  reasons  of  its  success  are  the  high  flexibility  of  this  method  to  adapt  to  different 
problems,  the  ability  to  deal  with  changes  of  topology  (contour  breaking  and  merging) 
without  any  extra  functions  and  the  guarantee  of  the  existence  of  solutions  in  the  class  of 
viscosity  partial  differential  equations  (PDEs).  Moreover,  extensive  numerical  algorithms 
based  on  Hamilton-Jacobi  equations  have  been  developed,  accurately  handling  shocks 
and  providing  stable  numerical  schemas. 

The  key  idea  of  the  LSM  is  to  implicitly  represent  a  contour  or  interface  as  the 
zero  level  set  of  a  higher  dimensional  function,  called  the  level  set  function  (LSF),  and 
formulate  the  evolution  of  the  contour  through  the  evolution  of  the  level  set  function.  For 
closed  contours,  signed  distance  functions  (SDFs)  were  originally  adopted  to  represent 
level  set  functions  because  they  directly  provide  stability  and  accuracy  to  the  LSM  [29,  30] . 
In  the  last  years,  however,  new  approaches  have  proposed  to  use  binary  functions,  rather 
than  distance  functions,  to  represent  the  contour.  This  change  of  representation  allows 
to  use  fast  and  convex  optimization  techniques  with  graph  cut  methods  [4,  3]  and  convex 
relaxation  methods  [7,  5,  16]  to  provide  powerful  alternatives  to  the  distance  preserving 
level  set  method.  Nevertheless,  distance  functions  are  still  essential  in  several  applications, 
e.g.  medical  image  segmentation,  surface  reconstruction  or  special  effects  in  computer 
graphics.  For  this  reason,  it  is  important  to  develop  a  fast  and  accurate  algorithm  for 
distance  preserving  level  set  methods. 

A  major  issue  of  distance  preserving  level  set  methods  is  the  limited  speed  of  the 
existing  algorithms,  which  are  based  on  Hamilton-Jacobi  equations  and  upwind  schemes 
[31,  30].  Two  main  issues  handicap  these  iterative  methods.  First,  the  speed  of  the  algo¬ 
rithms  is  limited  by  the  CFL  condition  [9],  which  is  a  necessary  condition  to  guarantee 
the  stability  of  PDE-based  iterative  schemes.  Secondly,  the  level  set  energy  is  hard  to 
optimize  because  of  the  f 1  -  based  total  variation  (TV)  term,  which  is  not  differentiable 
and  is  usually  regularized.  However,  the  regularization  significantly  slows  down  the  mini¬ 
mization  process  and  does  not  provide  an  exact  solution  to  the  problem.  In  our  work,  we 
will  overcome  these  speed  limitations  using  recent  efficient  techniques  in  E1  optimization. 

Another  major  issue  of  the  LSM,  pointed  out  by  Gomes  and  Faugeras  in  [18],  is 
a  contradiction  between  the  theory  and  the  implementation  when  LSF  are  represented 
by  SDF.  Indeed,  the  LSM  does  not  intrinsically  preserve  the  SDF  during  the  contour 
propagation  because  SDFs  are  not  solutions  of  the  Hamilton-Jacobi  equations  associated 
to  the  LSM.  Additional  techniques  are  then  necessary  to  preserve  the  LSF  as  an  SDF 
during  contour  evolution.  Unfortunately,  no  fully  satisfying  methods  have  been  proposed 
so  far.  The  two  most  common  approaches  that  have  been  suggested  to  fix  this  problem  are 
either  re-distancing  regularly  the  LSF  as  an  SDF  (this  procedure  is  called  re-distancing  or 
re- initialization) ,  or  constraining  the  LSF  to  remain  an  SDF  during  the  contour  evolution. 

Re-distancing  is  the  most  common  approach.  It  consists  in  stopping  the  evolution  of 
the  LSF  periodically  and  re- initializing  it  as  an  SDF  while  preserving  the  zero  level  set. 
This  approach  introduces  the  questions  of  when  and  how  to  re-initialize  the  LSF.  It  is  hard 
to  say  when  the  re-distancing  must  be  applied  as  there  is  a  trade-off  between  speed  (each 
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re-distancing  task  takes  time)  and  accuracy  (the  LSM  will  develop  irregularities  during 
the  evolution  without  re-distancing).  Similarly,  the  re-initialization  can  be  performed 
with  different  LSM  techniques,  using  PDEs  [36]  or  the  fast  marching  algorithm  [1,  29]. 
The  main  issue  with  the  re-distancing  approach  is  the  difficulty  to  preserve  exactly  the 
location  of  the  zero  level  set  during  the  re-distancing  process,  which  might  shift  the 
moving  interface  to  undesired  positions. 

The  second  approach  aims  at  constraining  the  LSF  to  stay  an  SDF  during  the  contour 
evolution,  avoiding  the  previous  re-distancing  procedure  altogether.  Gomes  and  Faugeras 
introduced  in  [18]  a  new  level  set  formulation  to  restrict  the  LSF  to  an  SDF.  The  new 
formulation  consists  of  three  coupled  PDEs,  which  makes  the  analysis  of  the  existence 
of  a  solution  and  the  numerical  implementation  more  difficult  than  the  standard  LSM. 
More  recently,  Li  et  al.  [24,  25]  proposed  to  add  a  penalty  term  in  the  level  set  energy  to 
constrain  the  LSF  to  be  close  to  an  SDF.  They  also  developed  a  new  algorithm,  simpler 
and  more  efficient  than  the  conventional  LSM  method.  However,  the  time  step  of  the 
algorithm  is  still  restricted  by  the  CFL  condition  and  the  SDF  property  is  only  encouraged 
but  not  enforced.  We  will  see  that  our  method  overcomes  all  of  these  limitations  and 
provides  an  efficient  way  to  constrain  exactly  the  LSF  to  be  an  SDF. 

The  level  set  method  was  first  introduced  in  computer  vision  to  carry  out  image  seg¬ 
mentation  [6,  21,  8]  and  then  extended  to  other  tasks  as  stereo  reconstruction  [11],  object 
tracking  [32],  and  object  recognition  [23].  In  this  paper,  we  focus  on  image  segmentation 
and  surface  reconstruction,  but  the  proposed  method  can  be  easily  extended  to  other 
problems.  In  image  segmentation,  the  level  set  method  re- formulates  the  parametric  ac¬ 
tive  contour  into  a  non-parametric  energy  minimization  problem  (i.e.  independent  of  the 
contour  discretization).  The  active  contour  is  embedded  as  the  zero  level  set  of  the  LSF, 
which  moves  according  to  the  active  contour  evolution  equation  and  drives  its  zero  level 
set  to  the  edges  of  the  desired  object.  The  level  set  formulation  easily  includes  edge-, 
region-  and  shape-based  terms  in  the  image  segmentation  criteria,  proving  the  high  flex¬ 
ibility  of  the  LSM  to  adapt  to  different  tasks  and  justifying  its  extensive  use  in  the  field. 
Similar  segmentation  models  have  been  applied  to  the  problem  of  reconstructing  a  surface 
form  unorganized  data  points  [41],  In  that  case,  the  LSM  allows  us  to  obtain  an  implicit 
representation  of  the  surface  by  an  SDF,  avoiding  complex  3D  parametrization  methods 
which  require  a  prior  knowledge  of  the  topology  of  the  surface,  and  providing  an  easy 
and  reliable  way  to  directly  estimate  surface  normals  and  curvature  from  the  level  set 
function. 

Here,  we  propose  a  fast  and  accurate  algorithm  for  distance  preserving  level  set  meth¬ 
ods.  We  constrain  the  LSF  to  remain  an  SDF  via  a  constrained  minimization  problem. 
The  constrained  minimization  problem  is  hard  to  solve  directly,  so  we  propose  to  split 
the  original  hard  problem  into  sub-optimization  problems  which  are  easier  to  solve,  and 
combine  them  together  using  an  augmented  Lagrangian  approach.  This  idea  is  borrowed 
from  the  split-Bregman  method  [17]  (and  more  generally  from  the  alternating  direction 
methods  of  multipliers  [12,  13,  26,  2,  35]),  which  is  an  efficient  El  minimization  method 
originally  developed  to  solve  the  compressed  sensing  problem.  The  proposed  algorithm 
holds  several  important  properties.  Firstly,  it  is  fast  because  it  is  not  limited  by  the 
CFL  condition,  leading  to  a  method  5  to  6  times  faster  than  most  distance  preserving 
LSM  algorithms  for  the  image  segmentation  problem.  This  improvement  is  due  to  the 
splitting  strategy  and  the  reformulation  of  a  constraint  minimization,  and  allows  us  to 
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deal  with  the  non- differentiability  of  the  TV  norm  and  go  beyond  the  CFL  time  step 
restriction.  Secondly,  our  algorithm  preserves  the  LSM  as  an  SDF,  avoiding  the  classi¬ 
cal  re-distancing  problem  and  providing  desirable  properties  for  some  applications.  For 
example,  this  makes  an  important  difference  in  surface  reconstruction,  where  surface 
normals  can  be  fast  and  reliably  estimated  during  the  surface  evolution  instead  of  be¬ 
ing  required  as  input  data  s.a.  [34,  22],  and  in  medical  image  segmentation,  where  the 
distance  information  can  be  exploited  to  include  topology  restrictions  into  the  problem. 
Finally,  our  algorithm  is  easy  to  implement  because  the  iterative  scheme  is  based  on 
standard  minimization  problems. 

The  rest  of  the  paper  is  organized  as  follows.  In  Section  2  we  review  very  briefly 
the  level  set  method  applied  to  image  segmentation  and  surface  reconstruction.  We 
introduce  our  algorithm  to  solve  efficiently  the  level  set  method  in  Sections  3  and  4. 
Section  5  presents  the  results  and  Section  6  draws  the  conclusions. 

2  Level  Set  Method  for  Image  Segmentation  and  Sur¬ 
face  Reconstruction 

Level  set  method  can  be  applied  to  perform  image  segmentation  using  the  active  contour 
method.  In  this  context,  the  segmentation  problem  is  defined  as  the  following  energy 
minimization  problem  w.r.t.  a  contour  C: 

min  f  Wb(s)ds  +  f  wlrn(x)dx+  f  w^ut(x)dx,  (1) 

CCn  JC  JCin  JCout 

where  0  is  the  image  domain,  s  is  the  arc-length  parametrization  and  the  first  term 
weights  the  length  of  C  by  an  edge  detector  function  uy,  such  as  in  [6,  21].  In  the  other 
terms,  C)n ,  Cout  designate  the  regions  inside  and  outside  the  contour  C  and  vj"1  ,  w°ut  are 
region-based  terms,  such  as  the  ones  proposed  in  [8].  We  adopt  this  minimization  model 
for  image  segmentation,  as  it  represents  a  large  class  of  active  contour  models  published 
in  the  literature. 

Actually,  similar  models  have  also  been  used  to  reconstruct  smooth  surfaces  from 
unorganized  data  points.  Given  a  set  of  (noisy)  points  {au}i<;<jv  lying  close  to  the 
unknown  surface,  a  smooth  estimate  of  the  surface  S  can  be  reconstructed  by  minimizing 
the  area  energy  weighted  by  the  distance  to  the  set  of  points  {ay}i<i<Ar,  i.e.  Wb(x)  = 
d{Xiy(x)  [20,  41].  Inspired  by  [22,  39],  we  propose  a  model  which  also  includes  a  region- 
based  term  to  improve  the  performance  with  sparse  data  sets.  We  adapt  the  model 
proposed  by  Lempitsky  and  Boykov  [22],  which  was  designed  to  align  the  normals  of 
the  reconstructed  surface  to  the  pre-computed  orientations  (n,:}i<,<Af  of  the  data  points. 
Based  on  this  surface  information,  a  semi-dense  vector  held  of  the  surface  normal  n  can 
be  estimated  everywhere  in  space  and  the  unknown  surface  can  thus  be  reconstructed  by 
iteratively  maximizing  the  alignment  of  the  reconstructed  surface  normal  N  ton.  Making 
use  of  the  divergence  theorem,  the  function  wlrn  can  be  identified  as  follows: 

argmax  /  n  ■  Afds  =  arg  min  /  —  di vndx,  (2) 

ccq  Jc  ccn  Jcin 


4 


which  corresponds  to  defining  wlrn  =  —  div  h.  The  flux  of  the  semi-dense  normal  field  at 
the  point  x  is  estimated  with  the  formula: 

El  I  sc  —  CCj|2 

e  2o-2  (x  —  Xi,  rii), 

i<i<N  V  2vrcr 

where  x— x*  denotes  the  vector  centred  at  x,  and  pointing  at  x  and  n,  is  the  surface  normal 
estimated  at  x%.  The  proposed  distance  preserving  LSM  can  estimate  the  surface  normal 
during  the  reconstruction  process:  rq  =  Vo  (ay),  where  0  is  the  SDF.  To  sum  up,  the 
proposed  functions  in  (1)  for  the  surface  reconstruction  problem  are:  wb{x)  =  d\Xi\{x), 
w™  =  —  div  ft,  and  w°ut  =  0.  Unlike  the  method  of  Lempitsky  and  Boykov,  the  proposed 
surface  reconstruction  algorithm  does  not  need  to  know  a  priori  the  normal  of  the  surface 
at  the  data  points. 

The  level  set  method  can  be  applied  to  solve  (1)  in  general  (and  specifically  the 
image  segmentation  problem  and  the  surface  reconstruction  problem),  by  re-writing  the 
minimization  in  terms  of  a  level  set  function  0  :  Q  — y  M  as  follows: 

min  [  wb(x)\V  H  (cf>)\  +  wr{x)H{4>)  s.t.  |V0|  =  1,  (4) 

*  Jn 

where  the  curve  in  M2  (or  surface  in  M3)  C  is  represented  by  the  zero  level  set  of  (f>,  H  is 
the  Heaviside  function,  and  the  constraint  V0|  =  1  guarantees  the  level  set  function  to 
be  a  signed  distance  function  [36]. 

3  An  Efficient  Algorithm  for  Level  Set  Method  pre¬ 
serving  Signed  Distance  Function 

In  this  section,  we  introduce  an  efficient  algorithm  to  solve  the  level  set  minimization 
problem  (4) .  The  main  idea  is  to  split  the  original  hard  problem  (4)  into  sub-optimization 
problems  which  are  well-known  and  easy  to  solve,  and  combine  them  together  using  an 
augmented  Lagrangian.  This  idea  is  borrowed  from  the  split-Bregman  method  [17],  which 
is  an  efficient  (Jl  optimization  method  recently  introduced  in  image  processing  to  solve 
the  compressed  sensing  problem. 

Let  us  consider  the  following  constrained  minimization  problem,  which  is  equivalent 
to  the  original  LSM  problem  (4): 


u=H((p) 

G 9=6 

q=vu  (5) 

p=V(f> 

IpI=i 

where  we  introduce  functions  <p(x),u(x)  G  M,  q(x),p(x)  €  Mn  (boldface  letters  are  used 
to  denote  vector  functions  and  n  =  2  for  2D  images,  and  n  =  3  for  3D  images).  The 
proposed  splitting  approach  makes  the  original  problem  (4)  easier  to  solve  because  (5) 
can  better  handle  the  non-differentiability  and  non-linearity  of  (4).  Indeed,  it  is  known 
from  [17,  38]  that  the  minimization  of  |V0|  can  be  carried  out  efficiently  by  decoupling 
the  £1-norm  |.|  and  the  gradient  operator  V.  The  term  |V0|  is  thus  replaced  by  |p|  and 


min  /  wb(x)\q\  +  wr(x)u  s.t. 

<j>,V,u,q,p  Jn 
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p  =  V0,  and  the  term  |Vi7(0)|  can  be  changed  into  |qf|,  <7  =  V«,m  =  H(p)  and  99  =  0, 
where  p  is  introduced  to  handle  the  non-linear  term  i7(0). 

Next,  we  want  to  reformulate  this  constrained  minimization  problem  as  an  uncon¬ 
strained  optimization  task.  This  can  be  done  with  an  augmented  Lagrangian  approach  [13], 
which  translates  the  constraints  into  pairs  of  Lagrangian  multiplier  and  penalty  terms. 

Let  us  define  the  augmented  Lagrangian  energy  associated  to  (5): 

C((f>,<p,u,q,p,  A)  =  /  wb\q\  +  wru 
Jn 

+  Ai  (>-0)  +  yO-0)2 
+  X2(u-H(p))  +  rj(u-H(p))2 
+  X3-(q-Vu)  +  7j\q-Vu\2 

+  A4-(p-V0)  +  ^|p- V0|2  s.t.  |p|  =  l  (6) 

where  A  =  (Ai,  A2,  A3,  A4)  are  the  Lagrangian  functions,  i.e.  Ai(x),  X2(x)  G  M,  A3(:r),  A4(.t)  G  Rn, 
and  ri, . . . ,  r4  are  positive  constants.  The  constraint  minimization  problem  (5)  reduces 
to  finding  the  saddle-point  of  the  augmented  Lagrangian  energy  £  [13].  The  solution  to 
the  saddle  point  problem  (6)  can  be  approximated  by  the  iterative  Algorithm  1. 


Algorithm  1  Augmented  Lagrangian  method  for  distance  preserving  level  set  methods 
1:  Initialize  0,  p,  u,  q,  p,  A 

2:  Find  a  minimizer  of  £  with  respect  to  variables  (cf),p,u,q,p)  with  fixed  Lagrange 
multipliers  Ak : 


(4>k,(Pk,uk,qk,pk)  =  argmin£  (cf),<p,u,q,p,Ak  x)  s.t.  |p|  =  1  (7) 

3:  Update  Lagrange  multipliers 

Ai  =Ai-1  +  ri((pk  -  4>k)  (8) 

\k=\k~l+r2{uk-H{pk))  (9) 

A3fc  =A3fc_1  +  r3(qk  —  S7uk)  (10) 

A/  =A/-1  +  r4(pk  -  V0fc)  (11) 


4:  Stop  the  iterative  process  when  ||0fc  —  0fc  1||2  <  e. 


We  initialize  0°,(/?O  with  the  signed  distance  function  of  the  initial  contour,  u°  = 
H(<p 0),  q°  =  V?/° ,  p°  =  V0°  and  the  Lagrange  multipliers  A0  =  0.  At  each  itera¬ 
tion,  an  alternating  minimization  method  is  used  to  find  an  approximate  minimizer  of 
£  (0,  (p,  u,  q,  p,  Afc_1)  considering  the  Lagrange  multipliers  Afc_1  fixed.  Then  the 
Lagrange  multipliers  are  updated  with  the  residual  associated  to  each  constraint  and 
the  process  is  repeated  until  the  change  of  the  level  set  function  0  falls  below  a  certain 
threshold  e,  which  happens  at  convergence. 
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In  general,  it  is  difficult  to  find  the  exact  minimizer  of  the  minimization  problem  (7) 
because  the  energy  (6)  is  not  convex  (there  is  no  guarantee  of  convergence).  However, 
experiments  show  that  a  good  approximation  can  be  found  by  the  alternating  direction 
method  of  multipliers  [13].  An  approximate  solution  is  thus  computed  by  iteratively 
alternating  the  minimization  of  C(<p,  u,  q.  p,  Afc_1)  w.r.t.  each  variable  while  considering 
the  others  fixed.  This  leads  to  Algorithm  2. 

The  next  step  is  to  determine  the  solutions  of  the  five  sub- minimization  problems  (12), (13), (14), (15) 
and  (16),  which  can  actually  be  computed  efficiently. 


Algorithm  2  Alternate  minimizations  for  an  approximate  solution  of  (7) 

1:  Initialize  0°  =  <p°  =  </A_1,  u°  =  uk~ 1,  q°  =  qk_1,  p°  =  pk~1. 

2:  For  l  =  1  and  and  fixed  Lagrange  multipliers  Ak,  solve  the  following  sub¬ 

problems  alternatively: 

ft  =argmin£(0,  <pl  1,v!  1,ql  1,pl  1,Al  x) 

4> 

(12) 

(pl  =  arg  min  £(0Z ,  <p,  ul~l ,  q 1-1 ,  pl~l ,  Al~l ) 

<p 

(13) 

ul  =  arg  min  £(</>*,  (p\  u,  ql~l,pl~l,  Al~l) 

U 

(14) 

ql  =  arg  min  C(4>1,  q,pl~l ,  Al~l) 

q 

(15) 

pl  =  argmin£(0z,  ul,  ql,p,  Az_1)  s.t.  |p|  =  1 

p 

(16) 

3:  Set  (cj)k,(pk, 

uk,qk,pk )  =  (4>L,  <pL,  uL,  qL,pL)- 

4  Sub-Minimization  Problems 

In  this  Section,  we  simplify  notation  by  omitting  the  super-index  and  the  tilde  symbol  in 
the  sub-minimization  problems  (12)-(16). 


4.1  Sub-Minimization  Problems  w.r.t.  0  and  u 

The  sub-minimization  problems  (12)  and  (14)  can  be  written  as  follows: 

T2 


mm 


/ 

Jn 


mm  / 

u  Jn 


T  2 

WrU  H - (  U  — 


rx)  ) 
^2 
r2)  ) 


—  7  y2  _|_  r_^_ 


V0-  [p 

Vu  —  (q 


The  Euler-Lagrange  equation  of  (17)  is: 

(— r4A  +  r4)0  =  — r4  divp 


Ai 

—  +  n  l  <p  +  — 


r\ 


r  i 


A4 

r4 

A3 

r3 


(V) 

(18) 


(19) 
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which  can  be  solved  efficiently  by  the  fast  Fourier  transform  (FFT)  as  proposed  in  [38,  37]. 
Let  us  denote  by  f(i,j )  a  scalar  2D  function  discretized  at  the  pixel  location  (i,j)  in 
the  discrete  image  domain  Q  =  [1,  Nx\  x  [1,  Ny],  We  also  define  the  identity  operator 
Xf(i,j)  =  f(i,j )  and  shifting  operators: 


Sxf(i,j)  =  /(*±  l,j) 

and  write  the  discretization  of  (19)  as: 

Syf(i,3 )  =  f(i,  3  ±  1)  > 

(20) 

[r4  («SX  -  2X  +  <S+  +  Sy 

whose  right  hand  side  is  discretized  as 

-  2X  +  S+)  -rj  <f>(i,j)  =  h(i,j ), 

(21) 

h(i,j)  =  -r4(X-Sx  )px{i,j )  -  r4(X-Sy  )py(i,j)+ 

~  (^  ~  S~)\4x(i,  j)  -  (I  -  Sy)X4y(i,  j)  +  j)  +  Ai(i,  j).  (22) 

Then  we  apply  the  discrete  Fourier  transform  J7,  taking  into  account  that  the  shifting 
operators  in  frequency  domain  correspond  to: 

977- 

FS\ tf  (yi,  Vi)  =  e±XZi7f  (yh  y5)  Zi  =  —yi  (23) 

2tt- 

FSff  (yi,  yj )  =  e±lZj 7 f  (yi,  yj)  Zj  =  —yj.  (24) 

where  (yi,Vj)  are  the  discrete  coordinates  in  the  frequency  domain.  Assuming  periodic 
boundary  conditions,  expression  (21)  in  Fourier  domain  now  equals: 

J2r4  (cos  (zj)  +  cos  (zj)  -  2)  -  nj  T<\>  (yi}  %•)  =  Th  (yi,  yj) ,  (25) 

Q 

and  provides  us  with  a  closed- form  solution  0*  of  (17): 

<f>*  =  F-1(Fh/g)  (26) 

which  we  can  efficiently  compute  by  FFT.  It  is  now  straightforward  to  apply  the  same 
procedure  to  compute  the  minimizer  u*  of  (18),  whose  Euler- Lagrange  equation  is: 

(— r3  A  +  r2)u  =  -wr  -  r3  div  q  +  — ^  +  r2  (^H(ip)  -  y-  j .  (27) 

4.2  Sub-Minimization  Problem  w.r.t.  (f 

The  sub-minimization  problem  (13)  can  be  written  as  follows: 

T  (t-v;))2 +ri(H^)~  (»+^))2  <28> 

Let  us  call  z  =  d>—  —  and  v  =  u  +  —  and  observe  that  the  minimization  is  decoupled  for 

r  r\  ri2  1 

each  pixel.  Let  us  define  the  function  we  want  to  minimize  as 


Observe  that  for  practical  implementations,  the  minimization  problem  (29)  involves  a 
smooth  approximation  H£  of  the  Heaviside  function.  We  propose  two  steps  to  find  quickly 
a  minimizer  of  (29). 

Step  1:  Find  a  solution  tp°  of  (29)  for  e  =  0  (i.e.  for  the  distributional/non- smooth 
Heaviside  function).  A  closed- form  solution  exists  for  this  problem  and  can  be  computed 
as  follows.  The  first  term  of  (29)  is  minimized  for  ip°  =  z.  As  the  distributional  Heaviside 
function  can  take  only  values  0  or  1,  the  second  term  of  (29)  is  minimized  for  ip°  <  0 
when  v  <  \  and  cp°  >  0  when  v  >  This  means  that  both  terms  are  minimized  for 
=  z  if  v  <  \  and  z  <  0  or  v  >  \  and  z  >  0.  Otherwise  we  must  choose  to  minimize 
the  greater  of  these  terms  and  set  ip°  =  0  if  E  (0)  <  E  ( z )  and  ip°  =  z  otherwise. 

Step  2:  Find  a  solution  ip°  of  (29)  for  e  >  0  using  the  standard  Newton’s  method  with 
ip°  as  initialization.  The  iterative  Newton’s  method  for  finding  the  minimizer  of  (29)  is 
as  follows: 

m+ 1  m  _  (0  m  _  O  (W™  ~  z)  +  r 2  ~  v)  4  (V™)  / 30^ 

^  ^  n  +  r2  -  v)  S's  ((pm)  +  r2Sj  (ipm)  ’  1 

With.  ipm=0  =  ip0, 

where  Se  and  8'e  are  smooth  approximations  of  the  Dirac  function  and  its  derivative. 
Finally,  we  have  observed  that  two  iterations  are  enough  to  find  a  good  approximation 
of  the  minimizer  of  (28),  leading  to  a  technique  almost  as  fast  as  a  closed- form  solution. 


4.3  Sub-Minimization  Problem  w.r.t.  q 

The  sub-minimization  problem  (15)  can  be  written  as  follows: 

min  f  wb\q\  +  ~  q  -  (Vu  -  — ) 

9  Jn  2  V  r3  / 


(31) 


Let  us  call  z  =  Vu  —  ^  and  observe  that  the  minimization  (31)  is  decoupled  for  each 
pixel.  The  solution  q *  of  (31)  is  given  by  soft-thresholding  with  the  shrinkage  operator 
[10]: 

q*  —  max  j|z|  —  — ,  o||~|  (32) 


4.4  Sub-Minimization  Problem  w.r.t.  p 

The  sub-minimization  problem  (16)  can  be  written  as  follows: 


min 

p 


\P\  =  1 


(33) 


Let  us  call  z  =  V0  —  ^  and  observe  again  that  the  minimization  (33)  is  decoupled  for 
each  pixel  and  is  equivalent  to  the  minimization  of  the  following  function  on  the  plane 
F(p)  =  f  |p  -  z|2  =  f  (|p|2  -  |z|2  -  2 p  •  z). 

We  can  introduce  the  constraint  in  |p|  =  1  in  F  by  writing  it  F(p)  =  y  (1  —  \z\2  —  2p  ■  z). 
It  is  obvious  that  the  minimizer  of  F  is  of  the  form  p*=  i.e.  when  vectors  z  and  p 
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have  the  same  orientation,  the  scalar  product  p  ■  z  reaches  a  maximum  and  the  constraint 
\p\  =  1  is  verified.  The  minimizer  p *  of  (33)  is  then  given  by: 

P  =  | r4V</>  -  A4| ~~  | r4V</>  -  A4|  A4 

5  Experiments  and  discussions 

Sections  5.1  and  5.2  present  the  results  in  image  segmentation  and  surface  reconstruction 
obtained  with  the  proposed  method.  Section  5.3  compares  our  algorithm  with  existing 
distance  preserving  LSMs. 

5.1  Image  segmentation  and  surface  reconstruction 

We  apply  Algorithm  2  for  image  segmentation  by  defining  the  edge  and  region  terms 
proposed  in  [6,  8],  which  have  been  extended  to  handle  both  gray  level  and  color  images. 
Figure  1  shows  the  results  obtained  for  different  images  from  the  Berkeley1,  Weizmann2 
and  GrabCut3  databases.  The  method  behaves  as  expected,  providing  the  same  results  as 
redistancing  or  penalty  methods  in  terms  of  final  segmentation,  but  with  a  considerable 
speed-up  in  time. 

We  also  use  the  proposed  model  to  successfully  reconstruct  several  surfaces  from  the 
Stanford  dataset4.  We  initialize  the  method  with  a  sphere  containing  all  the  data  points 
and  recompute  the  region  term  w™  =  —  div  h  every  5  iterations  as  described  in  Section 
2.  Unlike  [22,  39],  we  do  not  need  to  estimate  a  priori  the  surface  normal  at  the  data 
points,  as  the  surface  normal  at  the  points  is  directly  estimated  from  the  current  value  of 
the  LSF.  Finally,  the  process  is  sped  up  with  a  standard  multi-resolution  approach,  see 
Figure  2.  The  final  reconstructed  surfaces  are  shown  in  figures  3  and  4. 

5.2  Cortex  segmentation  with  coupled  surfaces 

Cortex  segmentation  [19]  is  a  problem  that  requires  the  distance  information  to  be  solved 
successfully.  Based  on  [18],  we  can  develop  a  segmentation  algorithm  for  the  cortical 
grey  matter  with  two  active  contours  coupled  by  their  relative  distance,  which  constrains 
the  thickness  of  the  cerebral  cortex.  Graph  cut  methods  and  convex  relaxation  methods 
cannot  be  directly  applied  to  solve  this  problem  because  they  use  the  binary  function  for 
representing  the  contour  and  binary  functions  do  not  hold  the  distance  information. 

The  cerebral  cortex  is  the  layer  of  the  brain  bounded  by  the  outer  and  inner  cortical 
surfaces,  that  is,  the  outer  interface  between  cerebral  spinal  fluid  (CSF)  and  grey  matter 
and  the  inner  interface  between  grey  and  white  matter.  Locating  this  cortical  surface  is 
a  first  step  in  many  brain  imaging  processes  and  measuring  its  thickness  is  a  common 
procedure  in  the  diagnosis  of  many  neurological  diseases.  We  will  see  that  the  use  of  SDF 
in  the  segmentation  of  cortical  surface  allows  us  to  include  information  about  the  cortical 

1http:  / /www. eecs.berkeley.edu/vision 

2http://www. wisdom. weizmann. ac.il/  vision 

3http://research. microsoft.com/en-us/um/cambridge/projects/visionimagevideoediting/segmentation/ 
grabcut.htm 

4http:  /  /  graphics.stanford.edu/data/3Dscanrep/ 
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Figure  1:  Proposed  level  set-based  segmentation  method  applied  to  natural  images.  For 
each  result,  we  show  the  segmentation  result  and  the  level  set  function.  We  plot  the 
initial  zero  level  set  in  blue,  the  final  contour  in  pink  and  the  ±1,±2  level  sets  of  the 
final  function  in  black. 
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( 


(a)  grid  30x30x30,  2.4s  (b)  grid  60x60x60,  10.6s  (c)  grid  120x120x120,  50.5s 

Figure  2:  Reconstructed  bunny  at  in  the  multiresolution  approach.  Linear  interpolation 
of  the  SDF  obtained  at  lower  resolutions  was  used  as  initial  LSF  for  higher  resolutions. 

structure  into  the  segmentation  problem  and,  at  the  same  time,  provides  an  estimate  of 
the  cortical  thickness. 

In  order  to  extract  the  cortical  layer  we  need  to  extract  its  two  bounding  surfaces  C 1 
and  C 2 .  In  theory,  therefore,  we  could  simply  detect  each  of  its  bounding  surfaces  or 
segment  the  regions  defined  by  the  white  matter  and  the  exterior  parts  of  the  brain  by 
independent  minimization  of  the  following  functionals  associated  to  C 1  and  C2 


(35) 


(36) 


where  Wb  is  a  boundary  indicator  and  vjv1  ,  uy2  are  region  descriptors  for  the  white  matter 
and  the  exterior  areas  of  the  brain.  In  practice,  however,  the  boundaries  between  grey 
and  white  matter  are  not  clear,  MRI  images  suffer  from  intensity  nonuniformity  and 
the  segmentations  obtained  with  local  region  descriptors  and  boundary  detectors  do  not 
correctly  locate  the  cortical  layer. 

To  overcome  these  limitations,  it  has  been  proposed  to  incorporate  a  constraint  on 
the  cortical  structure  to  obtain  a  better  segmentation.  For  simplicity,  we  modify  the 
model  proposed  in  [40]  to  suit  our  variational  formulation,  but  similar  models  for  cortex 
segmentation  s.a.  [28,  14]  could  also  be  adopted.  We  use  a  coupled  surface  model,  where 
a  functional  is  minimized  when  C 1  captures  the  CSF-grey  matter  interface,  C2  the  grey- 
white  matter  boundary  and  the  distance  between  them  is  close  to  the  expected  cortical 
thickness  d  (about  3mm).  The  problem  is  written  in  terms  of  LSM,  making  use  of  the 
SDF  0i  and  02  to  define  the  bounding  surfaces  C 1  and  C 1 .  The  functional  to  minimize 
is  then  given  by 


The  term  (0i  —  02  —  df  penalizes  segmentations  where  the  distance  between  the  bound¬ 
ing  surfaces  differs  from  the  expected  cortical  thickness  d.  Indeed,  when  0i  and  02  are 
the  SDF  defined  by  C 1  and  C2 ,  the  distance  between  the  surfaces  can  be  measured  on  the 
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(a)  grid  120x120x120,  50s 


(b)  grid  240x240x240,  391s 


(c)  grid  200x160x120,  151s 


(d)  grid  490x320x240,  1076s 


(e)  grid  200x160x120,  173s  (f)  grid  490x320x240,  1096s 


Figure  3:  Reconstructed  surfaces  from  scattered  data  points  at  different  resolutions. 
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(a)  grid  120x240x120,  173s  (b)  grid  240x480x240,  940s 


(c)  grid  60x60x120,  46s  (d)  grid  120x120x240,  315s 

Figure  4:  Reconstructed  surfaces  from  scattered  data  points  at  different  resolutions. 
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m 

(a)  Initialization,  (b)  Coupled  level  (c)  Non  coupled  (d)  Ground  truth,  (e)  k-mean  algo- 
set  functions.  level  set  functions.  rithm. 

Figure  5:  Segmentation  of  grey- white  matter  interface  on  MRI  images  of  human  cortex. 
The  segmentation  obtained  with  the  proposed  method  5(b)  is  clearly  closer  to  the  ground 
truth  5(d)  than  the  results  obtained  when  no  coupling  terms  is  considered  and  the  seg¬ 
mentation  is  performed  independently  for  the  inner  and  outer  cortical  surfaces  5(c).  The 
segmentation  5(e)  obtained  with  3  phases  of  the  k-mean  algorithm  fails  because  region 
descriptors  alone  cannot  segment  grey-matter. 


whole  domain  by  0\  —  02  and,  consequently,  the  term  {0\  —  02  —  d)2  drives  the  segmen¬ 
tation  to  solutions  where  the  distance  between  the  surfaces  is  consistent  with  the  cortical 
structure. 

The  minimization  technique  presented  in  Section  4  can  be  directly  applied  to  this 
problem.  The  same  splitting  variables  and  Lagrange  multipliers  are  now  defined  and 
solved  for  each  level  set  function  and  only  the  alternate  minimization  w.r.t  0i  and  02  are 
modified.  The  minimization  problems  w.r.t  0i  (the  analogous  applies  to  02)  is  now 


min 

<t>i 


,  U\2  ,  r4 

02~  d)  +  — 


(38) 


and  its  associated  Euler-Lagrange  equation 

(— r4 A  +  ri  +  c)0  =  —  r4  div p  +  —  +  rq  (ip  +  — )  +c(02  +  d)  (39) 

r4  \  rq  / 

can  be  solved  efficiently  by  the  fast  Fourier  transform  as  we  have  already  explained  in 
Section  4. 

We  have  applied  this  technique  in  an  illustrative  experiment  to  segment  the  cortical 
layer  in  different  slices  of  MRI  images.  Results  are  shown  in  Figure  5.  We  observe  that 
coupling  the  level  set  functions  and  constraining  the  expected  cortical  thickness  in  the 
segmentation  procedure  produces  a  segmentation  result  (Fig.  5(b))  close  to  the  ground 
truth  (Fig.  5(d)).  The  segmentation  obtained  without  coupling  (Fig.  5(c))  is  not  able  to 
find  the  fine  structures. 


5.3  Comparison  to  other  distance  preserving  level  set  methods 

In  this  section,  we  compare  our  proposed  algorithm  to  the  other  techniques  designed 
to  preserve  the  SDF  in  the  LSM.  To  that  purpose,  we  first  remind  the  necessity  of 
preserving  the  LSF  as  an  SDF  in  the  proposed  segmentation  task  with  a  simple  example. 
We  initialize  the  LSF  as  an  SDF  and  evolve  it  with  the  PDE  techniques  associated  to 
the  LSM  for  the  image  segmentation  problem  (4).  The  level  set  function  becomes  too 
steep  around  its  zero  level  set  after  a  few  iterations,  and  stops  its  evolution  because  the 
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(a)  Initial  contour  (b)  Initial  level  set  function,  (c)  Final  contour  de-  (d)  Final  level  set  function, 

defined  by  level  set  fined  by  level  set 

function.  function. 


Figure  6:  The  level  set  function  (0  =  0  in  pink  and  <p  =  ±1,±2  in  dark)  develops 
irregularities  during  the  propagation  and  stops  after  a  few  iterations  if  the  LSF  is  not 
constrained  to  be  (or  at  least  close  to)  an  SDF. 


geometric  features  (i.e.  curvature  and  normal)  are  not  correctly  estimated,  see  Figure  6. 
Two  common  approaches  have  been  introduced  to  overcome  this  problem,  either  by  re¬ 
distancing  the  LSF  or  by  maintaining  the  SDF  during  the  contour  evolution  (as  proposed 
in  our  method). 

We  will  compare  three  different  approaches  in  the  case  of  image  segmentation:  our 
method,  the  standard  re-distancing  approach  and  the  method  of  Li  et  al.  [24],  The  re¬ 
distancing  process  is  carried  out  with  the  fast  marching  method  [1] ,  while  the  method  of 
Li  et  al.  is  defined  by  introducing  a  penalty  term  in  the  energy  to  constrain  the  LSF  to 
be  close  to  an  SDF  as  follows: 

[  wb\S7H(<f>)\  +  wrH(<f>)  +  |(|V0|  -  l)2.  (40) 

We  refer  the  reader  to  [24,  25]  for  more  details. 

Our  algorithm  is  presented  on  Figures  7(c)-7(g).  Figures  7(h)-7(l)  show  the  re¬ 
distancing  method  [1],  Although  the  final  segmentation  and  the  final  LSF  provide  the 
desired  results,  the  periodic  re- initialization  process  produces  a  non- smooth  minimization 
of  the  level  set  energy  (observe  the  jumps  in  the  energy  plots).  Besides,  we  remind  that 
we  do  not  know  in  general  when  to  re- initialize  the  LSF  as  an  SDF.  In  our  experiments, 
we  applied  the  re-initialization  every  5  iterations.  We  also  remind  that  the  re-distancing 
process  is  not  guaranteed  to  exactly  preserve  the  location  of  the  zero  level  set  representing 
the  moving  interface. 

Next,  we  will  consider  the  method  of  Li  et  al.  [24],  which  is  more  related  to  our 
method.  Li  et  al.  introduced  the  nice  idea  to  constrain  the  LSF  to  be  close  to  an  SDF 
via  a  penalty  term.  However,  the  penalty  term  does  not  constrain  exactly  the  LSF  to  be 
an  SDF.  Besides,  the  value  of  the  penalty  constant  fi  is  a  trade-off  between  speed  and 
accuracy.  A  small  penalty  value  /i  in  (40)  provides  a  fast  algorithm  but  does  not  preserve 
faithfully  an  SDF,  leading  to  an  LSF  with  small  instabilities,  while  a  large  penalty  value 
//  provides  a  better  SDF  but  slows  down  significantly  the  minimization  process  as  the 
number  of  iterations  to  reach  the  convergence  state  increases  considerably.  Figures  7(m)- 
7(q)  present  the  results  with  a  small  value  /j  and  Figures  7(r)-7(v)  show  the  results  with 
a  large  value  /a. 
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Our  proposed  method  overcomes  the  limitations  of  Li  et  al.,  as  our  formulation  con¬ 
strains  the  LSF  to  be  an  SDF,  and  the  proposed  algorithm  is  fast  because  there  is  no 
need  to  assign  a  large  penalty  constant  to  the  penalty  term.  Indeed,  Figures  7(q),  7(v) 
and  7(g)  show  that  our  method  keeps  more  faithfully  the  LSF  as  an  SDF  because  the 
penalty  energy  is  lower,  by  at  least  one  order  of  magnitude  (the  minimum  of  the  penalty 
energy  is  around  1  with  our  method,  around  5  with  Li  et  aV s  method  with  a  large  /i,  and 
around  100  with  Li  et  a/.’s  method  with  a  small  //).  Our  method  is  almost  as  accurate  as 
redistancing,  when  it  comes  to  preserving  the  SDF,  and  overall  faster  than  redistancing 
and  Li’s  method,  as  shown  in  the  next  experiment  with  more  detail.  These  advantages 
are  provided  by  the  augmented  Lagrangian  approach,  which  can  preserve  accurately  the 
constraint  while  keeping  a  good  minimization  speed. 

To  quantify  the  improvement  obtained  with  our  method  with  respect  to  the  other 
distance  preserving  LSM  techniques,  we  have  used  the  previous  algorithms  to  segment 
72  images  from  the  Berkeley,  Weizmann  and  GrabCut  databases.  We  compare  on  Figure 
8  the  different  algorithms  in  terms  of  the  quality  of  preserving  the  distance  function  and 
speed.  Figure  8(a)  presents  the  values  of  the  penalty  terms  at  convergence  for  the  72 
segmented  images,  showing  that  our  method  preserves  the  SDF  almost  as  well  as  the 
redistancing  method  and  clearly  better  than  Li  et  a/.’s  method.  Indeed,  the  penalty 
values  for  our  method  and  the  redistancing  method  are  similar,  while  being  an  order  of 
magnitude  smaller  than  Li  et  a/.’s  method.  We  also  compare  the  time  for  each  method 
to  converge,  which  is  assumed  when  <  e.  We  observe  in  Figure  8(b)  and  Figure 

8(c)  that  our  algorithm  is  on  average  5  to  6  times  faster  than  the  redistancing’s  or  Li’s 
methods  for  all  images.  Finally,  we  present  on  Figure  8(d)  a  scatter  plot  of  the  penalty 
obtained  at  convergence  against  the  time  required  to  converge,  which  shows  that  our 
method  presents  a  good  trade-off  between  accuracy  and  speed. 

6  Conclusion 

We  have  introduced  an  efficient  algorithm  for  distance  preserving  level  set  methods  which 
overcomes  the  main  numerical  limitations  of  the  original  level  set  method,  i.e.  the  speed 
and  the  preservation  of  the  distance  function.  Although  other  fast  optimization  tech¬ 
niques  have  been  developed  for  the  level  set  method  [15,  7,  5,  16,  3],  they  cannot  preserve 
the  level  set  function  as  a  distance  function,  which  is  essential  in  some  applications  like 
surface  reconstruction,  medical  image  segmentation  or  segmentation  with  higher  order 
geometric  features  [27,  33].  Finally,  observe  that  the  proposed  algorithm  can  be  sped  up 
significantly  using  a  narrow-band  approach  and  a  parallelized  version  of  the  algorithm  on 
graphics  processing  units  (GPUs). 
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(1)  Redistancing  [1] : 
Penalty  energy 
w.r.t.  iterations. 


(q)  Li  et  al.  ’s  al¬ 
gorithm  [24]  (small 
penalty):  Penalty 

energy  w.r.t.  itera¬ 
tions. 


(r)  Li  et  al.  ’s  (s)  Li  et  al.  ’s  (t)  Li  et  al.  ’s  (u)  Li  et  al.  ’s  al- 

algorithm  [24]  algorithm  [24]  algorithm  [24]  gorithm  [24]  (large 

(large  penalty):  (large  penalty):  (large  penalty):  penalty):  Energy 

< f>  after  8, 000  </>  after  10, 200  Final  </>  after  evolution  w.r.t.  it- 

iterations,  8.98s.  iterations,  11.37s.  20,000  itera-  erations. 

tions,  22.02s. 


(v)  Li  et  al.  :s  al¬ 
gorithm  [24]  (large 
penalty):  Penalty 

energy  w.r.t.  itera¬ 
tions. 


Figure  7:  Level  set  method  with  our  algorithm,  with  the  re-distancing  procedure  [1],  and 
with  Li  et  aids  algorithm  [24].  The  first  three  columns  represent  the  evolution  of  the  LSF 
at  different  times  (the  zero  level  set  0  =  0  is^&i  magenta  and  the  iso  level  sets  0  =  ±1,  ±2 
are  in  dark).  The  fourth  column  plots  the  energy  fQ  Wb\ Vif(0)|  +wrH(cj))  +  |(|V0|  —  l)2 
and  the  last  column  shows  the  evolution  of  the  penalty  energy  J(|V0|  —  l)2  (closeness 
measure  between  LSF  and  SDF)  w.r.t.  iterations.  Our  method  provides  the  best  trade-off 
between  speed  and  preservation  of  the  distance  function. 
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Figure  8:  Comparison  of  quality  of  segmentation  and  speed  for  the  different  methods  on 
a  dataset  of  72  images.  Quality  of  the  LSF  is  measured  in  terms  of  the  penalty  term 
at  convergence  jV  f  (| VQ  —  l)2  when  the  obtained  contours  are  equivalent.  Our  method 
preserves  the  SDF  almost  as  well  as  redistancing,  clearly  better  than  Li  et  aV s  method 
and  is  faster  than  any  of  them. 
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